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Abstract 

The BogoHubov approximation is used to study the ground state and low- 
lying excited states of a dilute gas of atomic bosons held in an isotropic har- 
monic potential characterized by frequency u and oscillator length do- By as- 
sumption, the self-consistent condensate has a macroscopic occupation number 
A'o >> 1, with A^ — A'o << A'o- For negative scattering length — |a|, a simple 
variational trial function yields an estimate for the critical condensate number 

1 /2 

Nqc = (87r/25v^) ((io/|tt|) ~ 0.671 ((io/|tt|) at the onset of collapse. For pos- 
itive scattering length and large A'o >> do /a, the spherical condensate has a 
well-defined radius R >> do, and the low- lying excited states are compressional 
waves localized near the surface. The frequencies of the lowest radial modes 
(n = 0) for successive values of orbital angular momentum / form a rotational 
band uoi ~ ujqo + + 1) w (do/R)'^, with uoo somewhat larger than u. 

PACS numbers: 03.75.Fi, 05.30.Jp, 32.80.Pj, 67.90.+Z 
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Recent experimental demonstrations of Bose-Einstein condensation in dilute confined 
^^Rb [1] and ^Li [2] have stimulated theoretical research into their physical properties, 
based largely on the Bogoliubov approximation [3], originally introduced as a model for 
bulk superfluid ^He. Although this simple description of liquid He has long been familiar, 
much of its application to Rose condensed dilute atoms has involved extensive numerical 
analysis [4-6]. In contrast, Baym and Pethick [7] provide a more physical description of 
the confined ground state, emphasizing the relevant dimensionless parameters for ^^Rb, 
where the s-wave scattering length a is positive. The present work extends this picture 
to include negative values of a (as found, for example, in ^Li), along with the low-lying 
excited states for positive a. 

The Bogoliubov model is most simply understood by considering the familiar second- 
quantized field operators that obey boson commutation relations ['0(r), '(/;^^(r')] = 6{r — r'). 
The dynamics follows from the "grand-canonical hamiltonian" operator 

K = H-fxN = J dVi{;\T + U - fj,)i/; + 27rah'^m-^ J dVip^^i/;, (1) 

where H is the hamiltonian, N is the number operator, and /x is the chemical potential [8,9]. 
Here T = — 77.^V^/2m is the kinetic energy, U{r) is the external confining potential, and the 
short-range interatomic two-body potential has been approximated by a pseudopotential 
with an s-wave scattering length a. The presence of Rose condensation implies that the 
field operator has a macroscopic ensemble average = ^'(r), identified as the (in 

general, temperature-dependent) condensate wave function. For a dilute system at low 
temperature, most of the particles are in the condensate, and the deviation operator (f){r) = 
ip{r) — \I/(r) is treated as small (by definition, {4>{v)) — 0). An expansion of K through 
second order in these small field amplitudes immediately yields K « Kq + K', where 

Ko = J dV^*{T + U - fi)^ + 27rah'^m-^ J dVl^^l^, (2a) 
K' = y"cZy0t(T + C/-//)0 + 27ran^m-^y"dV(4|*|20V+*^0V^ + **^#), (2&) 



and the first-order contribution vanishes because the first variation of Kq provides the 
nonhnear Hartree equation for the condensate wave function [10,11] 

(T + t/ + 47ran^m-^|*|2* = 0. (3) 

In addition, the ensemble average of the total number operator N = Nq+N' determines the 
temperature-dependent number of particles in the condensate A^o — J dV and in the 
excited states N' = J dV ((^V)- The Bogoliubov approximation assumes that N' « Nq, 
thus neglecting terms of third and fourth order in the deviation operators; this assumption 
clearly fails sufficiently close to the onset temperature T^, since No(Tc) vanishes. 

The first step is to determine the condensate wave function ^, which then provides a 
static interaction potential for the low-lying excitations. Although the actual experimental 
traps are anisotropic [1,2], it is simplest to consider an isotropic three-dimensional harmonic 
potential U{r) = |ma;^r^, with a characteristic oscillator length do = ^JhJmu) (the effect 
of the anisotropy can be treated in perturbation theory). Following Baym and Pethick [7], 
I use a Gaussian trial function \l/(r) = C exp(— |r^/(i^), where C is a real normalization 
constant and the length scale d serves as the variational parameter. Substitution into 
Eq. (2a) yields the variational quantity 

Ko(/i, d) = ^7T^/^nuj[C'^{^ddl + fd^do"^) + C'^V27Tad^do^] - /iC^TT^/^d^. (4) 

The thermodynamic identity [8] dKo/dfj, = — A^o fixes the normalization = No/n^/'^d^, 
and the corresponding energy becomes 

E{X) = Ko + AtiVo = ^Nohuj[^{X^ + X'^) + aX^], (5) 

where X = do/d sets the spatial dimension of the spherical condensate, and 

a=^/2/^{Noa/do) (6) 

is a dimensionless parameter that characterizes the relative strength of the interparticle 
energy (note that a is proportional to the parameter = STrNoa/do = VS27r^ a defined 
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in Ref. [7]). In Eq. (5), the three terms represent the kinetic energy, the confining energy, 
and the interparticle energy, respectively. 

It is clear by inspection that E{X) becomes large for A ^ (large d/do) because of 
the spatial confinement in the harmonic potential, but the detailed behavior for large A 
(small d/do) depends on the value of the parameter a. In the absence of the interparticle 
interaction (u = 0), the minimum of Eq. (5) occurs at A = 1. For any positive a (repulsive 
scattering length with a > 0), the cubic term eventually dominates for A ^ oo, and the 
local minimum of E{X) remains absolutely stable for all cr > 0. For negative a (attractive 
scattering length with a < 0), however, the function E{X) diverges to — oo for A — > oo, and 
the local minimum disappears entirely at some critical negative value — jcrd, signaling the 
onset of an instability. 

The condition E'{Xo) = determines the location of the minimum, which satisfies 
the polynomial equation 1 = Aq + uAq. For |cr| << 1, the root is given approximately by 
Ao ~ 1 — |(J, with the corresponding energy £"0 ~ ^NohLv{l + ^a). As expected, a small 
repulsive (attractive) scattering length expands (contracts) the overall condensate size and 
raises (lowers) the overall energy. For large positive cr, it is straightforward to show that 
Ao ^ (J-^l^ - \ (7-^ with Eq ^ ^Nohuoa^/^, as found in Ref. [7]. 

The situation is very different for negative cr, because E{X) no longer has a global 
minimum, and even the local minimum disappears at the critical values Ac and ac deter- 
mined by the simultaneous conditions E'{Xc) = E"{Xc) = 0. An elementary calculation 
yields the values 

, 4 4 

Ac = 5^/^ 1.495, and ac = -— = — ^ ^ -0.535, (7) 

so that the interactions reduce the critical condensate size parameter dc ~ 0.669cio relative 
to that of the bare trap; the corresponding variational energy at the onset of the instability 
is Ec fti NqHuj. This calculation suggests that the energy becomes unbounded from 
below for a < ac through the disappearance of the local stable minimum rather than 
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through the onset of negative energy per particle. Since this variational calculation is 
merely an upper bound on the energy, the actual instability may well occur for less negative 
scattering lengths, and, indeed, numerical analysis from Ref. [5] gives the value Uc ~ —0.457 
for the vanishing of the ground-state energy per particle; the fa 15% difference in these 
values can be taken to characterize the accuracy of the variational estimate. 

As noted in Ref. [2] , this estimate predicts a maximum condensate number A^o ~ 1440 
for the parameters (a ~ —1.44 nm and do ~ 3.13 /xm) appropriate to the ^Li experiment. 
This value is an order of magnitude less than the total number of trapped ^Li atoms 
reported in Ref. [2]. Even at zero temperature, however, the nonzero interparticle poten- 
tial ensures that No{T = 0) < N, and finite-temperature excitation of quasiparticles (see 
below) further reduces No{T) below N; hence it is unclear whether this discrepancy repre- 
sents a failure of the Bogoliubov description. Although the effect of three-body clusters has 
recently been investigated [12], the small value of \a\ relative to the interparticle spacing 
suggests that two-body contributions dominate the physics of this many-body problem. 

The next step is to consider the noncondensate, which is described by the boson fields 
and 0t. Since K' in Eq. {2b) is a quadratic form in these field operators, it can be 
diagonalized by a canonical transformation [9], as in the original work of Bogoliubov [3] 
for a uniform condensate. Assume that the condensate wave function has the general form 
^'(r) = -\/]Vo e*'^^'"^ /(i")) where the real amplitude function / is normalized (J dV |/p = 1), 
and the phase S produces a particle current j = hNo\f\'^m~^'VS, as found, for example, 
in a singly quantized vortex [11,13,14], where the appropriate S is the angle in cylindrical 
polar coordinates. Define the linear transformation [9] 

0(r)=e^^« [ujir)a,-v*{r)al], 

' , (8) 

j 

where the primed sum means to omit the condensate mode from the sum. Here, the "quasi- 
particle" operators aj and ctj obey the usual boson commutation relations [q;^, ck^] = Sjk-, 

5 



[Q!j,Q:fe] = [q!],Q!^] = 0, ensuring that the transformation is canonical, and the wave func- 
tions Uj and Vj are chosen to satisfy the coupled "Bogoliubov" equations 

Luj — 4:7Ta?i^m~^\^\'^Vj = EjUj, 

(9) 

L*Vj — Anah m ^\^\'^Uj = —EjVj, 

where L — — (2m)~^/i^(V + iV/S)^ + [/ — /x + 87ra^^m~^|\E'|^ is a hermitian operator. It 
is easy to verify that the eigenvalues Ej are real and that the eigenfunctions obey the 
normalization J dV {u*Uk — v^Vk) = 5jk- Substitution of Eq. (8) into Eq. (2b) yields the 
simple and physical result [9] 

K' = -J2' E, I dV |^,f + E, a]a^, (10) 

so that the canonical transformation indeed diagonalizes the operator K'. In addition, if Uj 
and Vj are a solution with energy Ej , then the pair v* and u* are also a solution with energy 
—Ej-, since the quasiparticle number operator ajctj has nonnegative integral eigenvalues, it 
is necessary to take Ej > 0. Finally, Eq. (9) also has the solution uq = vq = f with Eq = 0, 
verifying that the Bose condensation does occur in the lowest self-consistent single-particle 
mode. Although these equations are easily rewritten in terms of two-component vectors 
(see, for example, pp. 477 and 501 of Ref. [8]), such formalism is unnecessary here. 

The structure of K' in Eq. (10) leads to a very simple description of the equilibrium 
states of the condensed Bose system. The quasiparticle ground state |0) satisfies the con- 
dition aj\0) = for all j ^ 0, and the excited states follow by applying arbitrary number 
of quasiparticle creation operators ct] to |0). In addition, the well-known properties of 
these harmonic-oscillator operators mean that the low-temperature properties are deter- 
mined entirely by the eigenvalues and eigenfunctions of the Bogoliubov equations (9). If 
(• • •) = Tr[- • • exp{—/3K') ]/Tr[ex.p{—/3K') ] denotes a self-consistent ensemble average at 
temperature T = {kBP)~^, then the only nonzero averages of one- or two-quasiparticle 
operators are (aja/e) = (cKfeaj) — 5jk = Sjkfj, where fj = [exp{PEj) — is the usual 
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Bose- Einstein distribution function. For example [9], the total number density n(r) has a 
condensate contribution no(r) = |\E'(r)p and a noncondensate contribution 

^'(r) = E' [f^ 1^^- Wl' + (1 + Z^-) 1^^- Wl'] ' (11) 

3 

where the condition N = No{T) + J dV n'{r) determines the temperature-dependent con- 
densate fraction No{T)/N. 

In the present case of a spherical condensate in a spherical confining potential U{r), 
where ^{r) = y/No f{f) satisfies Eq. (3), the Bogoliubov equations simplify greatly because 
the states can be characterized by the usual angular- momentum quantum numbers (l, m) 
associated with the spherical harmonics Yim-, along with a radial quantum number n. Given 
a solution for ^'(r), standard numerical techniques can determine the eigenvalues Enim and 
associated radial eigenfunctions Unim{f) ctnd Vnimii^) [6]- In order to gain more physical 
insight, however, it is valuable to consider a special limiting case in which the kinetic energy 
of the condensate wave function is negligible compared to the confining energy and the 
repulsive interparticle interaction energy. As discussed in [7] (see also Refs. [4] and [15]) 
this condition holds for a harmonic confining potential when the dimensionless parameter 
a = ■\/2fn [NQa/do) from Eq. (6) is large and positive, because the kinetic energy is then 
of order cr~^/^ relative to the other two contributions. As a result, the Hartree equation 
(3) for the condensate wave function then has the simple solution 

47ran^m-^|*(r)|^ = [ii - U{r)]^[^Ji - U{r)], (12) 

where 9{x) denotes the unit positive step function. If the oscillator length do is used to 
scale dimensionless lengths, the normalization condition on then yields the radius R of 
the spherical condensate 

= 15Noa/do = 15(7r/2)^/2 ^13) 

with chemical potential given by // = ^TkuR^. Although this approximation clearly fails in 
the immediate vicinity of the condensate surface (see, for example. Fig. 1 of Ref. [4]), its 
use in the Bogoliubov equations leads to only a small error in the limit u >> 1. 
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A combination of Eqs. (9) and (12) yields the following coupled eigenvalue equations 

[Dx + V{x)\ Unl{x) - V<{x) Vnl{x) = enlUnl, 

(14) 

-V^{x)Unl{x) + [Da; + V{x)] Vnl{x) = -e^lVnU 



where x = r /R and = 2Eni/huj. Here, 

(15a) 



1 d ^2 d + 



x"^ dx dx x"^ 

is the kinetic-energy operator, and V{x) = V<(a;) + Vy{x) is the potential energy, where 

V^(x) = R^{l-x^)d{l-x^), and Vy(x) = (x^ - 1) d{x^ - 1) (156) 

are both positive. Apart from the coupling between u and v, which occurs only for a; < 1 
through V<(a;), these equations look like those for radial eigenstates with orbital angular 
momentum I in an isotropic repulsive potential V{x) — i?^ |1 — a;^|, which has a central 
peak at the origin, vanishes linearly at a; = 1, and rises quadratically for a; >> 1. Thus 
the low-lying eigenfunctions are expected to be "surface" modes localized in the vicinity 
of the condensate surface at x = 1. 

In principle, these coupled differential equations can be solved numerically, but more 
physical insight comes from recognizing that they have a variational basis. If U{x) de- 
notes a two component vector with elements u{x) and v{x), then Eq. (14) has a matrix 
representation 

[Da: + Vix)]U{x) = ersUix), (16) 

where V^; = -Da; 1; V(a;) = V{x) 1 + V<(a;) ti, 1 is the 2x2 unit matrix, and are the 
familiar 2x2 Pauli matrices. It follows immediately that the variational quantity 

< J^x'dxUHx)[V, + V{x)]Uix) ^^^^ 
~ x^ dxW{x)T^U{x) 

provides an upper bound on the lowest eigenvalue e^i for each separate Z. As a very simple 

model, take 



with Jg°° dx \g{x)\'^ = 1. Substitution into Eq. (17) gives 



eo/ < A cosh 2x — B sinh 2%, 



(19) 



where 



A= x^ dxg{x)* [Dj; + V{x)]g{x) and 



S = 




Minimization with respect to x yields the condition tanh2x = B/A, with 



eoi < VA^ - 52. 



(21) 



If g also depends on some parameters, they can be varied to find the minimum of 
Eq. (21). For example, take g{x) oc x'^ 9{1 — x) + x~'^~^ 9{x — 1). The integrals in A 
and B are easily evaluated (with an integration by parts in the case of A), as is the 
normalization integral, and the minimum with respect to 7 was found numerically for 
several different values of I and R. Since the description is expected to hold best for 
larger R, only the case of i? = 5 will be considered in detail, corresponding to a value 
a ~ 168 for the dimensionless parameter a that determines the relative importance of the 
kinetic contribution in the energy balance for the condensate wave function. To a good 
approximation, the 11 lowest eigenvalues eoi for I = 0, • • • , 10 can be fit to a quadratic 
polynomial eoi ^ 5.83 + Z/25.05 + /^/25.25, which effectively has the intuitive form eoi « 
eoo + + ^) / R^ oi a radial zero-point energy eoo plus rotational energy + R"^ oi a rigid 
rotor; similar calculations for other values of R show that eoo(-R) depends only weakly on 
R, as expected from the form of the potential V{x). For comparison, the eigenvalues of the 
bare confining harmonic potential (here assumed isotropic) are 4n-|-2Z-|-3 in the same units 
of ^huj [16]. The most striking conclusion here is that the low-lying elementary excitations 
of the Bose condensate for relatively large condensate radius R and condensate number A^o 
should have a rotational band of states (those for n = and Z = 0, 1, 2, . . .) lying somewhat 
above the lowest state of the bare confining potential. With the values a ^ 10 nm and 
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do ~ 1.4 jim, which are appropriate for the experiment in Ref. [1], the radius R = 5 
corresponds to a condensate number A^o ~ 29, 000, about an order of magnitude larger 
than the value estimated in Ref. [1]. In Ref. [6], where the Bogoliubov equations were 
solved numerically in the presence of the "exact" condensate density, the corresponding 
dimensionless parameter is cr fa 9.7; thus it is not surprising that their eigenvalues (Ref. [6] 
reports only those for I < 4) differ considerably from the simple rigid-rotor form given 
above. 

The particle density operator p(r) = ip^ {r)tp{r) plays a central role in the response of a 
physical system to external perturbations. In the present case of a dilute Bose condensate, 
the noncondensate density p' = p — has an unusual and characteristic form 

p' ^ + (22) 

that follows from the Bogoliubov approximation introduced below Eq. (1). Since the oper- 
ators (f){r, t) and 0^(r, t) oscillate harmonically at the frequencies given by the eigenvalues 
of the Bogoliubov equations, Eq. (22) shows that the normal modes of the condensate 
can be identified as density (compressional) waves. In particular, the noncondensate part 
of the density-density correlation function becomes simply a correlation function of the 
deviation operators, given by 

thus a measurement of the frequency spectrum of density oscillations (for example, by 
studying the resonant response to small modulations of the trapping potential) would 
directly characterize the eigenvalues Ej. 
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